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Abstract: 

Quark-anti-quark {q-q) potentials with finite quark masses are studied from the q-q 
Nambu-Bethe-Salpeter (NBS) wave functions in quenched lattice QCD. With the use of a 
method which has been recently developed in the derivation of nuclear forces from lattice 
QCD, we derive the q-q potentials from the NBS wave functions. We calculate the q-q NBS 
wave functions in pseudo-scalar and vector channels for several quark masses. The derived 
potentials at each quark mass in both channels show linear plus Coulomb form. We also 
discuss the quark-mass and channel dependence of the q-q potentials. 
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1 Introduction 

An inter-quark potential is one of the most important ingredients of quantum chromo- 
dynamics (QCD). Experimentally, Regge slope [1] suggests that the inter-quark potentials 
show linear behavior at long distance. The string tension of the potentials between a quark 
and an anti-quark, a ~ l.SGeV, can be roughly estimated by hadron spectra using the re- 
lation J = M^/ {4a) with the spin J and the mass M of hadrons. At short distance, the 
inter-quark potential shows like the Coulomb interaction, which is, for example, suggested 
by the analogy between quarkonium and positronium. In fact, the linear plus Coulomb be- 
haviors of the inter-quark potentials reproduce the low-lying hadron spectra well in quark 
models. 

Theoretically, the study of the inter-quark potentials is challenging issue due to the 
non-perturbative nature of low energy phenomena in QCD. Lattice QCD simulation is the 
powerful tool for a numerical investigation in such a strong-coupling region of QCD. From 
the expectation value of Wilson loops, the potential for an infinitely heavy quark and anti- 
quark {Q-Q potential) and also the three-quark potential {3Q potential) can be obtained 
on lattices [2, 3]. The Q-Q potential from quenched lattice QCD simulations reveals the 
form of V{r) = ar - A/r with a = 0.89GeV/fm and A = 0.26. 

The actual inter-quark potentials suffer from the effect of quark motions, which is not 
included in the Q-Q potential. One can take into account the corrections coming from finite 
quark masses ruq order by order with the use of the heavy quark effective field theory. The 
effective field theory utilizes the hierarchy of scales coming from the heavy quark mass rUq 
and the relative velocity of heavy quarks, v. The potential nonrelativistic QCD (pNRQCD) 
is such an effective field theory at the ultrasoft scale uigv'^ obtained by integrating out the 
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hard scale rUg and the soft scale nigV [2, 4-7]. It is convenient to employ pNRQCD to 
obtain the corrections of heavy quarkonium spectra to the heavy quark motion. 

In this work, we study potentials between a quark and an anti-quark with a finite 
mass {q-q potentials) from quenched lattice QCD simulations. In order to explore the q-q 
potentials, we apply the systematic method which utilize the equal-time Nambu-Bethe- 
Salpeter (NBS) amplitudes to extract hadronic potentials [8-20] to the systems with a 
relatively light quark and an anti-quark. Due to the absence of the asymptotic fields of 
quarks, the reduction formula cannot be applied directly. Therefore, we assume that the 
equal-time NBS amplitudes for the q-q systems satisfy the Nambu-Bethe-Salpeter (NBS) 
equation with constant quark masses which could be considered as the constituent quark 
masses. By using the derivation of the relativistic three-dimensional formalism from the 
NBS equation developed by Levy, Klein and Macke (LKM formalism) [22-25], we shall 
obtain the q-q potentials without expansion in terms of nig. The preliminary results of the 
NBS wave functions and potentials of the q-q systems have been reported in Ref. [26]. 

The paper is organized as follows. In Sec. 2, we present our method to extract the 
q-q potentials. In Sec. 3, we show the lattice QCD setup. We then show our numerical 
results of the q-q wave functions and potentials in pseudo-scalar and vector channels for four 
different quark masses in Sec. 4. The obtained potentials reveal the linear plus Coulomb 
forms which are similar to the Q-Q potential from the Wilson loop. We perform fitting 
analyses of the q-q potential data. Sec. 5 is devoted to discussions and summary. 

2 Method of the extraction of inter-quark potentials 

Following the formulation to define the potentials on lattices [8, 9, 24], we show the basic 
equations to extract the q-q potentials on the lattice below. As shown in Ref. [24], the 
equal-time choice of the Nambu-Bethe-Salpeter (NBS) amplitudes satisfy the relativistic 
Schrodinger-type equation without an instantaneous approximation for original interaction 
kernels of the NBS equation. Therefore, we can start with the Schrodinger-type equation 
(which is reffered as LKM equation in Ref [24]) for the NBS wave function (/>(r) to define 
potentials: 



where fi{= mg/2) and E denote the reduced mass of the q-q system and the non-relativistic 
energy, respectively, and we simply assume nonrelativistic kinematics. Note that the po- 
tential U{f,r') is generally energy-independent and non-local [9, 21]. In appendix A, we 
discuss potentials derived from the relativistic Schrodinger-type equation. The relativistic 
effects may be necessary to reproduce the meson mass spectra in the wide energy region. 

For the two-nucleon case, it is proved that the Schrodinger-type equation is derived 
by using the reduction formula [9]. Due to the absence of asymptotic fields for confined 
quarks, we suppose that the q-q systems satisfy the NBS equation with their constant 
quark masses. In this study, constant quark masses rUg are determined by half of vector 
meson masses My, i.e., rUg = My/2, as usually taken in constituent quark models. Then, 
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one finds Schrodinger-type equation of Eq. (2.1) as a three-dimensional reduction of NBS 
equation by applying LKM [22-25] method. 

The energy-independent and non-local potential U{f, r') can be expanded in powers of 
the relative velocity v = —iV/pi of q-q systems at low energies, 

U{r,r') = V{r,v)5{r-r') 

= {Vlo (r) + Vnlo (r) + ---)6{r- ?) , (2.2) 

with 

VLo{r) = Vc{r) + VT{r)Si2. (2.3) 
VNLo{r) = VLs{r)L-S, (2.4) 

where the N^LO term is of order 0{v^), and 5i2, L and S being the tensor operator, 
orbital angular momentum and spin of the q-q systems, respectively. Note that the velocity 
expantion is different from the usual l/niq expansion, and the central force Vc{r) in the 
leading order includes not only linear and Coulomb terms but also higher order terms in 
the l/mq expansion such as the spin-spin interaction, Vsp\a{f)(yq ■ Sq, which is regarded 
as an order 0{l/m'^). This spin-spin interaction is an important ingredient of the mass 
formula in the constituent quark model. At the leading order, one finds 

V{r) VLoir) = ^ + E. (2.5) 

The effective leading order q-q potentials Vlo(^) a.re studied in this work. The convergence 
of the expansion of v can be checked by studying the energy dependence of the local poten- 
tial Vlo{^ as in Ref. [12]. When the local potential Vlo{^) has little energy dependence, 
the potential between q-q is well described only by Vlo(^)- In contrast, if the energy de- 
pendence is large, higher order terms are necessary. The study of the energy dependence 
of the local potential in the q-q systems is an important future work. 

In order to obtain the NBS wave functions of the q-q systems on the lattice, let us 
consider the following equal-time NBS amplitudes 

x{x + f,x,t-to;J'') 
= {0\q{x + f,t)rq{x,t)Jqq{to;J^) |0) 

= ^An{0\q{x + r,t)rq{x,t)\n) e~*^"(*-*»), (2.6) 

n 

with the matrix elements 

An = {n\Jqq{to;r)\0). (2.7) 

Here T represents the Dirac 7-matrices, and Jqq{tQ\ J"") stands for the source term which 
creates the q-q systems with spin-parity on the lattice. The NBS amplitudes in Eq. 
(2.6) are dominated by the lowest mass state of mesons with the mass Mq at large time 
separation (t ^> to): 

x{r,t-tQ-,r) = -^^x(^ + ^=',^,i-to; J'') 

^ Ao(/>(f; J^)e-*^«(*-*o), (2.8) 
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with V being the volume of the box. Thus, the q-q NBS wave function is defined by the 
spatial correlation of the NBS amplitudes. 

The NBS wave functions in S-wave states are obtained under the projection onto zero 
angular momentum [P^''^^^), 

^(^■^ •^'^) = ^ E r), (2.9) 

geo 

where g E O represents 24 elements of the cubic rotational group, and the summation is 
taken for all these elements. Using Eq. (2.5) and Eq. (2.9), we will find the q-q potentials 
and NBS wave functions from lattice QCD. 



/3 


a lattice size 


volume 


-^conf 


6.0 


0.104 fm 32^ X 48 


(3.3fm)3 


100 



Table 1. Simulation parameters used in this work. Scale is set by string tension [3, 27]. 
3 Numerical setup of the lattice simulations 

In this section, we show the actions and simulation parameters in this work 
the standard plaquette gauge action, 

X,fJ,,U 

with /3 = 2Nc/g'^. The plaquette P^y is defined as 

P^^u = U^{x)U,{x + fi)Ul{x + y)Ul{x), (3.2) 

where U^{x) is a link variable. As for quark fields ip{x), we adopt the standard Wilson 
fermion action, 

SfI^P, i^,U]^Yl i^(^)K{x, vMy), (3.3) 

x,y 

K{x,y) = 

Sx,y - K ^{{1 - J^,)U^ix)5cc+il,y + (1 + 7A.)f^^(y)4,j/+/i}, (3.4) 

where k is the hopping parameter. 

We generate the quenched gauge fields on a 32'^ x 48 lattice with QCD coupling 
/3 = 6.0, which corresponds to the physical volume V = (3.3fm)^ and the lattice spac- 
ing a = 0.104fm [3, 27]. We measure the q-q NBS wave functions for four different hopping 
parameters k = 0.1520, 0.1480, 0.1420, 0.1320: the corresponding pseudo-scalar (PS) me- 
son masses Mps in the calculation are 0.94, 1.27, 1.77, 2.53GeV, and vector (V) meson 
masses Mv=1.04, 1.35, 1.81, 2.55GeV, respectively. The number of configurations used 



We employ 
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Figure 1. The q-q NBS wave functions in PS(a) and V(b) channels. The wave functions are 
normahzed at origin. All the wave functions are localized in the box and indicate the bound states. 

in this simulation is 100 for each quark mass. The simulation parameters are summarized 
in Table 1. The calculation of the q-q NBS wave functions requires gauge fixing, because 
q and q operators are spatially separated at the sink time slice. Here we adopt Coulomb 
gauge, which is frequently used for studies of hadron spectroscopy in lattice QCD. As for 
the source operator of the q-q systems, we employ a static wall source in Eq. (2.6), 

Jggito;Jn=Qito)TQ{to), (3.5) 

with the static wall quark operator 

Qito)=Y,^(^^to). (3.6) 

X 

We note that gauge fixings and sink operators can be arbitrary chosen in the formalism, 
and all potentials with different gauge fixings and sink operators give the same physical 
observables, i.e., mass spectra and scattering lengths for instance. In this work, we employ 
Coulomb gauge and a local operator for sink operators. We can take another gauge fixing 
and sink operator, and the potential obtained with these conditions is generally different 
from that obtained in this work. In appendix B, we discuss the sink-operator dependence 
of the potential by using a gauge-invariant smeared operator for the sink. 

4 Numerical results for the q-q potentials 

First, we show the numerical results of the NBS wave functions in Fig. 1. Fig. 1(a) and 
(b) are the NBS wave functions for each quark mass in PS and V channels, respectively, 
at the time slice t = 20. The NBS wave functions mostly vanish at r = 1.5fm for all quark 
masses in both channels. This indicates that the spatial volume V = (3.3fm)^ is large 
enough for the present calculations. The size of a wave function with a lighter quark mass 
becomes larger than that with a heavier one. Comparing the results in PS and V channels, 
little channel dependence between PS and V channels is found, although the quark-mass 
dependence of the wave functions is a bit larger for V channel. 
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Figure 3. Plots of the potential with arbitrary constatnt energy shift V{r) — E = V^(/)(r)/(2/i(/)(r)) 
in PS channel (a) and V channel (b) for each quark mass. The solid curves are the fit function with 
linear plus Coulomb form shown in Table 2. 



In Fig. 2, we show the Laplacian parts of q-q potentials in Eq. (2.5), V^(/)(r)/(/)(r), for 
each quark mass and channel. Fig. 2(a) shows V'^ cl){r) / (j){r) = 2fi(V{r) — E) in PS channel 
for each quark mass at the time slice t = 20. As shown in Fig. 2, one can see that the 
potential form is similar to that obtained from Wilson loop, namely, the potential form 
looks like linear plus Coulomb form, although the derivation of the potentials is largely 
different between these two methods. Fig. 2(b) represents V^(^(r)/(/>(r) in V channel for 
each quark mass at the same time slice i = 20. The basic properties are similar to that in 
PS channel, although quark mass dependence is a bit larger for V channel. 

Fig. 3(a) and (b) show the potentials with arbitrary energy shifts E, i.e, V{r) — 
-E=V^(/)(r)/(2//(/>(r)) in PS and V channels, respectively, for each quark mass at the time 
slice t = 20. Note that the quark mass mq(= 2/i) is determined by the half of vector meson 
mass, ruq = My/2, as mentioned in the previous section. 

We perform fit analyses of the potentials in Fig. 3(a) and (b). For the fit function, we 
choose the linear plus Coulomb form, /(r) = ar — A/r+C. We fit /(r) to the potential data 
for each quark mass and channel. We use the on-axis data with the range 3 < r/a < 10 in 
the fit. The fit results are summarized in Table 2, and denoted by solid curves in Fig. 3(a) 
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and (b). x^/-^df is around 0.5 for all the fit, which means the data are well described by the 
linear plus Coulomb form. We find moderate quark mass dependence of the string tension. 
The string tension becomes larger as increasing quark masses in both channels, and that for 
heaviest quark mass in our simulation is about 820MeV/fm, which is comparable to that 
obtained from an expectation value of the Wilson loop. On the other hand, the Coulomb 
coefficient A strongly depends on quark masses. The Coulomb coefficient becomes small as 
increasing quark masses, and is roughly approaching to that obtained from an expectation 
value of the Wilson loop. 







Pseudo-scalar 






Vector 




K 


cr{mg,i) 


Aimg,PS) 




a{mg,i) 


A{mg,V) 






MeV/fm 


MeV-fm 




MeV/fm 


MeV-fm 




0.1320 


819(47) 


215(7) 


0.32 


825(48) 


195(7) 


0.63 


0.1420 


753(34) 


264(5) 


0.35 


765(37) 


216(6) 


0.61 


0.1480 


691(30) 


338(5) 


0.46 


723(39) 


249(7) 


0.44 


0.1520 


601(29) 


443(5) 


0.31 


697(63) 


291(13) 


0.23 



Table 2. The fitting results of the potentials in Fig. 3. The function to be fitted is f{r) = 
a{mq,i)r — A{mq, i)/r + C{mq,i). The fit range is 3 < r/a < 10. 

Next, we perform another type of fit analyses. Assuming that the string tension 
cr is independent of the quark masses due to the quenched QCD simulations, where the 
contributions from quark loops are eliminated, we perform the fit by minimizing the general 
xV^df which is defined as xV^df = E™,,* X^("i<7, 0/^df with i = PS, V [30]. We caU 
the fit "universal fit" here. The fit functions /i(r) can be explicitly written by /i(r) = 
ar — A{mq, i) /r+C{mq,i). The free parameters of the fit are fi, A{mq, i), C{mq, i) for fi{r). 
In the fit, we choose the range of the potential data as 3 < r/a < 10. The fitting results 
are shown in Table 3. The general x^/-^df is achieved with 1.52 for /i(r) with a = 723 (30) 
MeV/fm. Since our simulation includes all the quark mass effects, /i(r) is modified by the 
higher order effect of l/niq expansion. 





Pseudo-scalar 


Vector 


K 


A{mq,PS) 


Aimq,V) 




MeV-fm 


MeV-fm 


0.1320 


231.45(4.50) 


212.42(4.60) 


0.1420 


270.63(4.53) 


224.91(4.64) 


0.1480 


331.91(4.95) 


248.96(5.00) 


0.1520 


414.93(5.89) 


284.76(7.03) 



Table 3. The fitting results of the potentials for /i(r) without the channel and the quark mass 
dependences of the string tension. The obtained string tension is cr = 723(30) McV/fm. The fit 
range is 3 < r/a < 10. 
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We also perform another type of the universal fit with the fit function, /2(?") = ar — 
A{mg,i)/r + C{mg,i) + Bi{mq,i)log{r) — B2{mg,i)/r'^ , which contains l/rUg correction 
terms, i.e., and log(r). The log term is predicted by pNRQCD and the effective string 
theory, respectively [28, 29]. In this fit, we use not only on-axis data but also off-axis 
data, and results are summarized in Table 4 (x^/-^df = 2.14). As shown in Table 4, if we 
employ the fit function f2{'f) in which 0(1 /niq) terms are taken into account, the Coulomb 
coefficients are smaller than those obtained fi{r) and are comparable to the values from 
Wilson loop. 







Pseudo-scalar 






Vector 




K 


A{mg,PS) 
MeV-fm 


Biimg,PS) 
MeV 


B2img,PS) 
MeV-fm2 


A{mg,V) 
MeV-fm 


Biimg,V) 
MeV 


B2img,V) 
MeV-fm2 


0.1320 
0.1420 
0.1480 
0.1520 


93.6(139) 
90.5(172) 
93.2(142) 
91.0(65.4) 


161(128) 
164(91.9) 
123(76.2) 
48.2(112) 


20.8(43.4) 
26.0(42.7) 
45.2(35.8) 
78.8(19.7) 


103(129) 
81.3(123) 
101(106) 
276(267) 


119(128) 
173(97.8) 
146(97.2) 
1.66(44.2) 


21.2(38.0) 
15.9(34.8) 
19.2(35.5) 
1.13(78.2) 



Table 4. The fitting results of the potentials for the fitting function /2(?') = or — A{mq,i)/r + 
C{mq,i) + Bi{mq,i)\og{r) — B2{mq,i)/r^ . The obtained string tension is (t = 705(120) MeV/fm. 
The fit range are is 5 < r/a < 12. 



5 Discussion and summary 

We have studied the inter-quark potentials between a quark and an anti-quark {q-q po- 
tentials) from the q-q Nambu-Bethe-Salpeter (NBS) wave functions. For this purpose, we 
have utilized the method which has been recently developed in the calculation of nuclear 
force from QCD [8, 9]. We have calculated the NBS wave functions for the q-q systems 
with four different quark masses in pseudo-scalar and vector channels and obtained the 
q-q potentials through the Schodinger-type equation. In this framework, the q-q potentials 
basically contains full quark motions with the finite masses. As a result, we have found 
that the shapes of the q-q potentials are the linear plus Coulomb form which is similar to 
the static Q-Q potential obtained from Wilson loop. 

For the fitting, we have employed two types of fitting functions. One is the linear plus 
Coulomb form regarded as the leading order (LO) terms in the l/rUg expansion. The other 
function includes the next leading order (NLO) terms in addition to LO terms. The fitting 
results with LO terms reveal that the Coulomb coefficients depend on the quark masses and 
are larger than those predicted from Wilson loop. On the other hand, if we have employed 
the NLO terms together with the LO terms, the Coulomb coefficients become smaller and 
are comparable to the value from Wilson loop. With the both fitting functions, we have 
obtained the string tension which is comparable to the value from Wilson loop. 

This is the first step to study the q-q potentials from the NBS wave functions, and 
the main purpose of the present study is to show that the method is applicable to the q-q 
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Figure 4. Comparison of "relativistic" potentials V^' (r) witli non-relativistic potentials l/^-f^- (r) 
shown in the main part. 



potentials. We find that the obtained q-q potential has the basic property of that obtained 
from Wilson loop. Therefore, this method can be used for the study of the q-q potentials 
with finite quark masses. 
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A Relativistic effect for the potential 

In the main part of the paper, we assume the non-relativistic kinematics for the Schrodinger- 
type equation. We here calculate the potential with relativistic kinematics to evaluate the 
relativistic effect for the q-q systems. 
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The relativistic Schrodinger-type equation for the q - q system in the continuum limit 
is written as 

+Vir}(Pif) = E(P{f), (A.l) 

with the leading order potential V{'r) of the velocity expansion. On a lattice, a discrete 
Fourier transformation of r gives sin(p). Thus yjp''^ + nT-q is replaced by ^ {sm{j?))'^ + m^, 
and the integral becomes a summation on lattice. Similar as the procedure in the main 
part, we obtain the potential from the discretized version of the relativistic Schrodinger- 
type equation of Eq. (A.l). 

Fig. (A) is a comparison of potentials with relativistic kinematics, V^'{r) and non- 
relativistic potential y'^-^-(r) shown in the main part. The relativistic potentials for k = 
0.1320 (k = 0.1520) are shown by open square (filled square) in Fig. (A), while the non- 
relativistic potentials for k = 0.1320 (k = 0.1520) are shown by open circle (filled circle). 
As shown in Fig. (A), V^'{r) show the linear plus Coulomb behavior, and V^-{r) and 
V^'^'{r) is similar for large r. On the other hand, for small r, the difference between them 
becomes large as naturally expected, and the short distance differences of the potentials 
may contribute to reproducing higher energy states. 

We calcurate the root mean square radius (r^) from the NBS wave function for 
the heaviest q - q system (k = 0.1320) employed in this work. Then we obtain (r^) = 
0.3550(14) fm, and the relativistic potential in Fig. (A) is consistent with non-relativistic 
one around 0.36 fm. Thus, for the heavy quark case such as the charm quark, we conclude 
that the non-relativistic potential picture gives the proper description of the heavy quarko- 
nium. On the other hand, for the light quark such as strange quark, short range part of 
the relativistic potential differs from non-relativistic one. This shows that the relativistic 
effects starts to contribute to q - q systems. 

B Potential from a gauge-invariant smeared operator 

As we showed, the obtained potentials exhibit a Coulomb plus linear behavior. However, 
potentials with different operators are generally different. Therefore, the Cornell-like be- 
havior is not universal. Here, we show a potential with a different operator from that used 
in the main part. 

The NBS amplitude with gauge-invariant smeared sink operators is defined by 

= (0| q{x + f, t)L{r, x, t; m)rq{x, t)Jgg{to; r) |0) , 
L{r = nfi, x, t) 

= Uf,{x + nfL,t) ■ ■ ■ Uf,{x + fi,t)Uf,{x,t) (B.l) 

A schematic figure of the amplitude is shown in Fig. (5). The operator L{f,x,t;m) con- 
structed by link variables connects x and i?-|-r with a straight-line path. Here, the direction 
of r is chosen to be the x-, y-, or z-direction, i.e., on-axes. 
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Figure 5. Schematic figure of gauge-invariant smeared operator. 



Figure (6) shows a potential obtained from the smeared NBS amphtude, V^'^^{r) (red 
points for V channel and blue points for PS channel), and that obtained in Coulomb gauge, 
yConi.^^^ (green points for V channel and blue points for PS channel). Note that the data 
of V^^^{r) are only calculated on the points with integral multiples of the lattice spacing a, 
because r in Eq. (B.l) is on-axis. V^^^{r) shows the linear plus Coulomb behavior similar 
to that in Coulomb gauge, and, more over, the two potentials almost coincide. This fact 
shows that the gauge-invariant operator is also a suitable one for a constituent quark mass, 
and the Coulomb-gauge operator used in the main part is similar to the gauge-invariant 
operator of Eq. (B.l). 




Figure 6. Comparison of a potential with gauge-invariant smeared operator, F™'^(r), to that 
with Coulomb-gauge operator, y'-^°"'-(r). The red (blue) points are the data with gauge -invariant 
smeared operator in V (PS) channel, and the green (magenta) points are that in Coulomb gauge in 
V (PS) channel. The two potentials almost coincides. 
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